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1.  Introduction 


Regional  seismic  monitoring  and  discrimination  capabilities  that  are  desirable  under  a 
potential  Comprehensive  Test  Ban  Treaty  (CTBT)  can  be  improved  by  developing  new 
algorithms  and  procedures  for  distinguishing  between  earthquakes,  nuclear  explosions  and 
mining  explosions  of  various  kinds.  Much  effort  in  past  discrimination  studies  has  concen¬ 
trated  on  extracting  various  features  of  the  spectrum  that  are  characteristic  of  earthquakes, 
nuclear  explosions  or  mine  blasts. 

One  particular  spectral  feature  that  characterizes  some  mining  explosions  is  a  modu¬ 
lation  of  the  spectrum  introduced  by  a  ripple-fired  explosion.  A  ripple-fired  event  usually 
involves  detonation  of  a  number  of  explosions  that  are  often  regularly  grouped  in  space  and 
time.  Such  explosions,  known  as  quarry  blasts,  have  low  magnitudes  that  may  be  close  to 
those  of  nuclear  explosions  that  one  might  monitor  under  the  CTBT.  As  examples  of  these 
kinds  of  mine  blasts,  we  consider  using  array  data  from  the  Arctic  Experimental  Seismic 
Station  (ARCESS)  in  northern  Norway,  previously  analyzed  by  Der  et  al  (1993).  The  top 
panel  of  Figure  1  shows  a  single  channel  from  a  typical  mine  blast,  sampled  at  40  points 
per  second,  with  four  arrival  phases  identified.  The  arrival  phases  shown  correspond  to 
body  waves ,  which  are  subdivided  into  compressional  P-waves,  denoted  here  by  Pn  and  Pg, 
and  shear  S-waves ,  denoted  by  Sn  and  Lg.  At  regional  distances,  Pn  is  usually  the  first 
to  arrive,  having  traveled  down  through  the  crust  and  then  horizontally  near  the  top  of 
the  upper  mantle.  Pg  is  a  wave  that  comes  in  later  than  Pn,  and  travels  the  path  between 
source  and  receiver  wholly  within  the  Earth’s  crust.  Sn  is  a  shear  wave  that  travels  a 
path  similar  to  that  of  Pn  but  arrives  later  because  shear  waves  propagate  slower  than  the 
P-waves  that  make  up  Pn.  Lg  is  a  type  of  guided  shear  wave  that  has  most  of  its  energy 
trapped  in  the  Earth’s  crust.  In  the  bottom  panel,  we  show  the  first  arrival  phase  at  all 
five  channels. 

The  pattern  induced  by  the  delay-firing  that  is  characteristic  of  mining  activity  can 
be  represented  in  terms  of  a  model  that  contains  a  signal  s(t )  and  its  echoes,  observed  at 
a  set  of  time  delays,  say  t*,  T2, . . . ,  rn,  producing  a  theoretical  output  of  the  form 

n 

x(t)  =  s(t)  +  J2  eosit  -  rj)>  ( 1 ) 

3= 1 

where  the  initial  signal  s(t)  appears  at  the  set  of  time  delays  indicated,  with  amplitudes 
0i,  02,  •  •  • )  #n-  Identifying  the  length  of  the  time  delays  for  arrival  phases  as  well  as  their 
configuration  will  be  helpful  for  distinguishing  ripple-fired  patterns  from  those  due  to  path 
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effects  or  local  reverberations.  A  number  of  authors  have  examined  various  aspects  of 
this  problem  and  have  proposed  techniques  for  analyzing  these  ripple-fired  seismic  signals. 
Chapman  et  al  (1992)  show  reflection  patterns  for  a  number  of  delay-fired  configurations 
and  propose  a  cepstral  deconvolution  method  for  estimating  the  delays.  Baumgardt  and 
Ziegler  (1988)  consider  lining  up  the  log-spectra  and  cepstra  for  an  array  and  looking 
for  common  reflection  patterns.  Alexander  et  al  (1995)  extend  this  analysis  by  adding 
up  or  stacking  the  cepstra.  Hedlin  et  al  (1990)  propose  graphical  techniques  involving 
threshold  modifications  of  the  time  varying  log-spectra  and  cepstra.  During  the  course  of 
this  project,  we  have  developed  three  different  approaches  based  on  two  different  signal 
and  noise  models.  They  are: 

(i)  deconvolution  of  a  common  signal  on  an  array  (see  Shumway  and  Der,  1985,  Der 
et  al,  1987, 1992,  1993),  assuming  an  additive  noise  model,  followed  by  a  cepstral 
analysis  of  the  deconvolved  signal, 

(ii)  a  search  through  Box-Jenkins  seasonal  autoregressive  moving  average  ( ARMA) 
models  with  parameters  equal  to  the  duration  and  time  delays  in  a  multiplicative 
signal  and  noise  model 

and 

(iii)  a  nonlinear  cepstral  F-statistic  (see  Shumway  et  al,  1997)  based  on  a  frequency 
domain  multiplicative  model  for  the  signals  and  noises.  The  approach  formalizes 
ideas  of  Alexander  et  al  (1995)  and  Baumgardt  and  Ziegler  (1988)  that  are  based 
on  the  premise  that  a  common  reflection  pattern  should  appear  at  each  channel 
on  the  array. 

The  deconvolution  approach  in  (i)  uses  an  additive  linear  model  for  the  observed  series 
at  the  jth  sensor  in  an  array,  say  yj(t),  j  =  1, 2, . . . ,  N  that  takes  the  form 

yj(t)  =  a,j(t )  <8>  x(t )  +  nj(t),  (2) 

where  the  notation  <g>  denotes  convolution  and  the  impulse-response  functions  aj(t)  are  used 
to  denote  the  path  and  receiver  effects  that  distinguish  the  jth  element  of  the  array  from 
others.  It  will  be  assumed  that  these  functions  have  relatively  smooth  Fourier  transforms  so 
that  they  can  be  estimated  from  the  spectra  of  the  observed  series.  We  assume  throughout 
that  the  signal  and  noise  are  stationary  and  independent  of  each  other.  Given  reasonable 
estimators  for  the  path  and  receiver  transforms,  the  multiple  deconvolution  procedure  of 
Shumway  and  Der  (1985)  is  used  to  develop  the  deconvolution  estimator  for  the  common 
function  x{t).  Then,  because  of  (1),  we  may  use  the  cepstral  analysis  procedure  of  Bogert, 
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Healy  and  Tukey  (1962),  denoted  by  BHT  in  the  sequel,  to  estimate  the  time  delays  for 
the  ripple  firing.  Section  2  gives  details. 

It  is  also  possible  to  consider  a  multiplicative  noise  model  of  the  form 

Vj ( t )  =  a,j (t)  ®  x(t)  <S>  rij ( t )  (3) 

as  an  alternate  representation,  leading  to  procedures  (ii)  and  (iii)  mentioned  above.  A 
simplification  of  (3)  with  equal  receiver  effects  and  delays  Tj  =  jd,j  =  1, 2, . . . ,  n  propor¬ 
tional  to  some  underlying  delay  d  allows  recasting  (3)  as  a  seasonal  autoregressive  model 
of  the  form  considered  in  the  classic  work  of  Box  and  Jenkins  (see  Box,  Jenkins  and  Rein- 
sel,  1994  for  expository  material).  The  delays  are  the  seasonal  periods  and  the  duration 
or  number  of  reflections  corresponds  to  the  order  of  the  seasonal  moving  average  term. 
A  low-order  autoregressive  component  models  the  combined  effects  of  source,  path  and 
instrument  response.  Seasonal  ARMA  models  are  searched  over  a  number  of  plausible 
delays  and  duration,  with  the  best  value  of  a  Bayesian  information  criterion  (BIC)  due  to 
Schwarz  (1978)  used  to  select  the  best  model.  This  approach  is  covered  in  Section  3. 

Finally,  the  Fourier  decomposition  of  the  observed  spectrum  in  the  multiplicative 
model  (3)  suggests  Method  (iii)  above.  In  this  approach,  the  observed  spectrum  is  ex¬ 
pressed  as  a  product  of  signal  and  noise  spectra  and  Fourier  transforms.  Taking  loga¬ 
rithms  yields  an  additive  model  in  the  logarithms  of  the  squared  Fourier  transforms  that 
we  can  use  as  a  model  for  detecting  a  common  set  of  periodic  components  on  an  array  of 
suitably  detrended  log-spectra.  In  our  approach,  assuming  a  multiplicative  noise  model, 
detrended  log-spectra  are  derived  as  realizations  of  stationary  processes  whose  periodic 
signal  components  are  frequencies,  with  periods  proportional  to  delay  time  differences. 
Using  an  approach  proposed  by  Shumway  (1971)  for  detecting  a  common  signal  in  a  col¬ 
lection  of  stationarily  correlated  series,  an  F-Statistic  is  derived  that  is  proportional  to  the 
stacked  cepstrum  and  the  spectrum  of  the  stacked  log-spectra.  Advantages  of  this  cepstral 
F-Statistic  stem  from  its  superior  resolving  power  and  the  fact  that  statistical  significance 
can  now  be  asserted  for  selected  delay  peaks. 

Simulated  array  data  and  data  involving  four  phases  extracted  from  ten  Kola  mining 
explosions,  measured  at  ARCESS  in  Scandinavia,  are  used  to  compare  the  time  and  fre¬ 
quency  domain  approaches.  We  cover  the  three  methodologies  in  Sections  2,  3  and  4  and 
the  data  analysis  in  Section  5. 
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2.  Multiple  Deconvolution 


The  primary  aim  of  the  multiple  deconvolution  procedure  is  to  isolate  a  reasonable  estima¬ 
tor  for  the  unobserved  stationary  process  x(t)  in  the  additive  model  (2).  This  problem  is 
somewhat  different,  than  single  channel  deconvolution,  where  there  is  only  one  realization 
for  the  signal  x(t).  In  what  follows,  we  suppose  that  the  spectra  of  the  signal  and  noise 
processes  are  denoted  by  Ps{y )  and  Pn{v)  respectively  with  v,  -1/2  <  v  <  1/2  defined 
as  frequency  in  cycles  per  point.  We  assume  the  signal  to  noise  ratio  to  be  known,  so 
there  are  no  parameters  to  be  estimated  by  maximum  likelihood  (see,  for  example,  Han¬ 
nan  and  Thomson,  1974  and  Shumway,  1988  for  approaches  that  develop  parametric  and 
non-parametric  estimators  for  the  signal  and  noise  spectra).  In  the  present  case,  it  is  pos¬ 
sible  to  show  (see  Shumway,  1988)  that  the  optimal  deconvolution  estimator  for  x(t)  is  of 
the  form 

N 

x(t)  ~  *52  h3  (*)  ®  Vj (*)  (4) 

j  = 1 


where  hj{t)  are  filters  applied  at  each  sensor,  with  frequency  response  functions  determined 


as 


H,{v)  = 


A-Jv) 


Ejlil  + 


(5) 


where  82{u)  =  Pn(v) /Px(v)  and  Aj(v)  is  the  Fourier  transform  of  aj(t).  Since  the  signal 
and  noise  spectra  are  somewhat  complicated,  as  shown  in  the  next  paragraph,  we  use  a 
small  ridge  correction  52(u)  =  82,  assumed  to  be  constant  over  frequency,  with  the  constant 
taking  the  role  of  the  usual  regularization  parameter  used  to  stabilize  solutions  to  poorly 
conditioned  inverse  problems.  The  resulting  estimator,  x(t),  should  consist  of  the  common 
components  of  all  array  elements  and  should  be  free  of  the  effects  of  the  path  and  receiver. 
Note  that  we  still  must  know  values  for  the  frequency  response  functions  Aj  (v) ;  we  discuss 
this  problem  later. 


Now,  from  (1),  the  spectrum  of  x(t)  is  given  by 


PxH  =  I  B(v)\2Pa(v), 


(6) 


with 

n  n 

\0{y)  |2  =  X!E  6i6k  C0S[27rI'(rj  ~  Tk )],  (7) 

j= 0  fc=0 

( 0o  =  1)  which  is  the  component  introduced  by  the  possible  ripple  firing.  BHT(1962)  noted 
that  the  squared  transfer  function  (7)  will  be  periodic,  with  components  proportional  to 
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the  time  delay  differences.  For  multiplicative  models  like  (6),  BHT  proposed  an  additive 
model,  obtained  by  transforming  to  logarithms.  For  a  sample  series  x(t),  t  =  0, 1, . . . ,  T- 1 
and  discrete  frequencies  of  the  form  vt  =  £/T,t  =  0,1,...,  T/  2,  the  log  spectrum  can  be 
written  as 

log  Px(vi)  =  log|0(z/£)|2  +  log  P8(vt).  (8) 

BHT  noted  that  the  first  component  will  be  periodic,  with  frequencies  proportional  to 
T/ijj—Tk),  j,  k  =  1, . . . ,  n  and  proposed  computing  the  spectrum  of  (8),  called  the  cepstrum 
which  should  exhibit  the  periodicities  in  (8)  as  cepstral  peaks. 

The  results  given  above  suggest  deconvolving  on  the  array  to  get  the  estimator  for  x(t) 
and  then  computing  the  cepstrum  of  the  deconvolved  series  to  isolate  the  periodicities.  As 
mentioned  before,  the  success  of  the  procedure  will  depend  on  finding  reasonable  estimators 
for  the  receiver  response  functions  Aj( v).  To  do  this,  we  note  that  the  spectrum  of  the 
received  signal  in  (3)  will  be 


Py,  M  =  l4fMI2-P.M  +  Pn{v).  (9) 

We  suppose,  furthermore,  that  the  noise  spectrum  can  be  neglected  in  (9)  and  that  the 
function  Aj(v)  is  reasonably  smooth  over  frequency.  Consider  the  cubic  spline  approxima¬ 
tion 

log  \Aj{vi)\2  +  log Px(v)  «  aj o  4-  anvt  +  aj2 v}  +  aj3 vj  +  ajA{vt  -  vf)%,  (10) 

where  uf  is  the  knot  location  and  {y-vf)\  for  l  =  0, 1, . . . ,  T/2.  The  sphne  approximation 
with  one  knot  should  be  smooth  enough  so  that  the  residuals  from  the  fitted  spline  will 
still  retain  the  periodic  components  in  \6{v)\2  that  correspond  to  the  ripple  fired  signal. 
The  resulting  approximation  (10)  for  the  smooth  part  of  log  |  Ay(z/)  |2  can  be  exponentiated 
and  the  positive  square  root  of  the  result  used  as  the  frequency  response  function  Aj(u). 

To  illustrate  some  of  the  above  ideas,  consider  the  simulated  ensemble  of  signals  shown 
in  the  top  panel  of  Figure  2.  We  generated  different  second-order  autoregressive  processes 
with  a  frequency  content  comparable  to  the  seismic  waveform  in  Figure  1  and  modulated 
them  with  a  fixed  exponentially  decaying  function  that  emulates  a  transient  signal.  The 
results  were  added  together  at  time  delays  r0  =  0,  T\  =  8,  t2  =  15,  r3  =  23,  r4  =  31  with 
amplitudes  9$  —  1, 6\  =  .9, 62  =  .9, 03  =  .6,64  =  .7  to  simulate  a  typical  realization, 
according  to  (1).  The  five  resulting  series,  shown  in  Figure  2,  all  contain  the  slightly 
different  signals,  at  the  same  set  of  delays  and  we  can  argue  that  the  model  given  by  (1) 
and  (2)  is  plausible  and  that  any  method  based  on  the  model  should  work  reasonably  well. 
A  typical  observed  spectrum  is  shown  in  the  bottom  panel  of  Figure  2  and  we  note  that 
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it  appears  to  be  the  sum  of  a  smooth  trend  component,  log  \Aj(u)\2,  fitted  as  a  cubic  spline 
with  a  single  knot  at  0.25  cycles  per  point,  and  a  periodic  component,  corresponding  to 
log|0(z/)|2. 

Figure  3  shows  the  bank  of  deconvolution  filters  hj(t),j  =  1,...,5  in  the  left  panel 
and  the  resulting  deconvolved  signal  x(t)  in  (4).  To  find  the  cepstrum  of  the  deconvolved 
signal,  take  the  discrete  Fourier  transform 

T— 1 

X(i)  =  T"1/2  53  x(t)  exp{-27T ilt/T]  (11) 

£= 0 

of  the  deconvolved  series  and  compute  the  spectrum  of  the  log  periodogram  as 

T/ 2-1 

C(d)  =  T~1/2  53  log|i'(^)|2exp{27 rUd/T},  (12) 

£=0 

expressed  as  a  function  of  delay,  d,  in  points.  The  cepstrum  of  the  deconvolution,  shown  in 
the  bottom  right  panel  of  Figure  3,  isolates  delays  at  r  =  7, 8, 15, 23,  and  29  points.  From 
the  nature  of  the  contrived  data,  we  expect  delays  at  7, 8, 15, 16, 23,  and  31  points  with 
potentially  the  strongest  components  corresponding  to  time  delays  of  7  or  8,  appearing 
four  times,  15  or  16,  which  appear  three  times,  23  points  which  appears  twice  and  31 
points  which  appears  once. 

3.  Seasonal  Multiplicative  ARMA  Models  for  Arrays 

There  is  a  form  of  (3)  that  fits  into  the  multiplicative  seasonal  ARMA  developed  by 
Box  and  Jenkins  (see  Box,  Jenkins  and  Reinsel,  1994).  For  this  case  we  let  x(t)  in  (1) 
be  defined  for  time  delays  rk  =  kd,  k  =  1, 2, . . . ,  n  that  are  multiples  of  some  underlying 
assumed  firing  delay  d,  so  that  one  might  envision  a  multiplicative  model,  driven  by  white 
noise  processes  Wj  (t) ,  of  the  form 

Vj(t)  =  ^l^Wj(t)  (13) 

for  the  output,  where  B  denotes  the  delay  operator,  i.e.  Bkx(t )  =  x(t  -  k),  so  that  the  n 
reflections  at  multiples  of  d  are  expressed  by  the  polynomial  operator 

n 

e(Bd)  =  Y^0kBkd  (14) 

k-Q 
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(15) 


with  $o  —  1.  The  autoregressive  operator  of  order  p, 

k=l 

defines  the  input  signal.  Dargahi-Noubary  (1995)  has  argued  that  many  seismic  source 
function  models  can  be  regarded  as  having  been  generated  by  the  one-parameter  third- 
order  autoregressive  model  (p=3),  with 

<t>{B)  =  {l-e-koBf  (16) 

where  k0  is  an  unknown  parameter.  Combining  (13)  and  (16)  leads  to  a  seasonal  ARMA 
with  a  third  order  autoregressive  component  and  a  dth  order  moving  average  component 
with  period  n  that  can  be  written  as 


Vj  W  ~  3e  *°  yj  (t  -  1)  +  3e  2fc°  Vj  (t  -  2)  -  e~3k° Vj  (t  -  3)  =  Wj  ( t )  +  hwj  ( t  -  kd) .  (17) 

1 

In  the  above  model,  the  delay  d  corresponds  to  the  spacing  of  the  charges  in  the  ripple  fired 
event  and  the  number  of  pulses  n  might  be  proportional  to  the  number  of  shots.  Hence,  the 
model  selected  will  depend  on  both  d  and  n,  whereas  the  parameters  estimated  under  each 
model  are  k0, 6U  02, . . . ,  6n  and  the  variance  of  Wj(t),  say  a^,  assuming  a  common  variance 
for  the  stochastic  inputs  generating  the  signal.  The  underlying  time  domain  model  is  more 
specific  in  specifying  the  duration  or  seasonal  increment  and  in  assuming  that  the  smooth 
component  of  the  spectrum  is  generated  by  a  single-parameter  third-order  autoregressive 
model.  Note  that  the  general  third-order  model  has  also  been  proposed  by  Tjpstheim 
(1975)  as  a  model  for  short  period  seismic  data. 


Minimizing  the  sum  of  squared  errors  for  specified  values  of  n  and  d  leads  to  the 
objective  function 

N  T 

SSE(k0,9i,e2,...,6n )  = E  E  «#),  (is) 

j=l  t=nd 

which  maximizes  the  log  likelihood,  and  we  may  estimate  a2  as 


=  i v^EE^w-  (is) 

3  t 

where  the  nonlinear  optimization  only  involves  NT'  —  N(T  -  nd+1)  residuals.  In  order 
to  select  a  model,  we  choose  n  and  d  as  the  joint  minimizers  of  the  Bayesian  Information 
Criterion  (BIC)  of  Schwarz,  1978,  say 


BIC(n ,  d)  =  logo-2  + 


(n  + 1)  log(NT') 
NT' 


(20) 
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The  Schwarz  criterion  is  consistent,  i.e.,  the  probability  of  fitting  the  wrong  order  model 
converges  to  zero,  whereas  procedures  based  on  Akaike’s  Information  Theory  Criterion 
(AIC)  (Akaike,  1973),  are  not  and  lead  to  overfitting  in  some  cases.  If  one  has  the  idea 
that  there  exists  a  true  model  of  some  order  and  the  sample  is  quite  large,  as  in  the  present 
situation,  it  is  more  reasonable  to  use  the  Schwarz  criterion.  If  one  is  inclined  to  think 
of  finding  the  best  approximation  to  some  unknown  model  of  possibly  infinite  order,  the 
various  AIC  measures  are  asymptotically  efficient  (Shibata,  1980)  whereas  BIC  is  not.  For 
a  summary  of  the  nonlinear  Gauss-Newton  estimation  procedure  applied  to  the  repeated 
measures  ARMA  model,  see  Shumway(1988). 

To  give  an  example,  consider  searching  the  contrived  data  in  Figure  2  for  the  best 
fitting  seasonal  ARMA  process.  It  is  convenient  to  limit  the  number  of  reflections  to  the 
possible  range  1  <  n  <  6  and  the  delays  to  the  range  4  <  d  <  13.  We  begin  with  d  =  4,  so  as 
not  to  confuse  the  delays  with  the  first  three  lags  of  the  autoregressive  part.  Figure  4  shows 
the  resulting  values  of  BIC  and  we  note  that  there  are  a  number  of  local  minima,  mostly 
occurring  at  d  =  8.  The  global  minimum,  however,  is  attained  for  n  =  4,  d  =  4,  and  would 
lead  to  values  that  differ  from  the  correct  values  n  =  4,  d  =  8.  The  estimates  6j,j  =  1,4  in 
this  latter  case  were  -.47,  .58,  -.60,  .16  and  the  alternating  signs  are  not  reasonable  for  the 
application.  The  most  reasonable  model,  for  n  =  3,  d  =  8  had  the  third  smallest  BIC  and 
gave  positive  estimates  (.83,  .41,  .15)  for  the  parameters  and  k0  =  .98.  The  ambiguity  in  the 
multiple  local  minima,  largely  due  to  the  difficulty  of  estimating  the  number  of  reflections 
n,  suggests  looking  at  values  of  BIC  averaged  over  n,  shown  in  the  lower  panel  of  Figure 
4.  The  consistency  of  the  minima  for  d  =  8  suggests  taking  this  value,  which  turns  out  to 
be  reasonable  from  the  individual  time  delays,  placed  at  T\  =  8,  t-j  =  15, 7-3  =  23,  T4  =  31 
points  respectively  with  weights  6X  =  .9,  02  =  .9, 63  =  .6, 64  =  .7.  In  summary,  there  are 
more  possible  interpretations  for  the  ultimate  model  implied  by  the  time  domain  analysis 
and  few  approximate  statistical  significance  tests  are  available  for  the  number  of  reflections. 
One  might  consider  testing  against  a  model  with  n  =  0,  but  we  did  not  do  so  because  of 
the  number  of  possible  alternatives. 
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Stacked  Plot  of  BIC  for  Simulated  Data 


Stacked  Plot  of  Average  BIC  for  Simulated  Data 


Figure  4:  BIC  Contours  for  Seasonal  ARMA  Search  Using  Contrived  Data.  The 
darkest  value  corresponds  to  the  minimum.  The  bottom  panel  averages  BIC  over  the 
number  of  reflections  and  shows  a  consistent  delay  at  d  =  8  points. 
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4.  The  Nonlinear  Cepstral  F  Detector 

We  may  also  take  a  nonlinear  approach  based  on  the  multiplicative  model  (2).  Note  that 
the  spectrum  of  the  observed  series  yj(t)  in  this  case  has  the  representation 

P„»  =  |A 

=  |J4i(r)|2|«(I/)|2PsMP„JM> 

using  (6).  Then,  taking  logarithms  yields  the  additive  decomposition 

log Pyj(u)  =  log  |-4j(^) |2  +  log  \9{y)\2  +  log Pa(v)  +  log Pnj(v) 

To  follow  up  on  the  suggestion  in  Section  2,  based  on  the  additive  decomposition  above, 
we  consider  computing  the  logarithm  of  the  spectral  estimator  at  each  channel.  Defining 
the  discrete  Fourier  transform  of  yj(t)  as  Yj  (£)  (see  (11))  at  frequencies  v£  =  £/T,£  = 
0, 1, . . . ,  T/2  —  1,  we  have  the  representation 

log  \Yj(£)\2  =  log  Ps{vt)  +  log \Aj{ve)\2  +  log  \9(u£)\2  +  log  \Nj(£)\2,  (21) 

where  \Nj(£)\2  is  the  periodogram,  known  to  have  approximately  a  chi-squared  distribution 
with  2  degrees  of  freedom  when  the  noise  rij(t),  for  example,  is  a  linear  process  (see 
Hannan,  1971,  p.  249).  Subject  to  regularity  conditions  on  the  linear  process,  the  variance 
is  proportional  to  to  the  square  of  the  spectrum,  with  an  error  term  that  is  0(T-1). 
Furthermore,  for  two  frequencies  separated  by  multiples  of  1/T  ,  the  correlation  is  of  the 
same  order.  Since  the  variance  is  proportional  to  the  square  of  the  mean,  we  may  expect  the 
logarithm  to  have  an  approximately  constant  variance.  If  one  smoothes  the  periodogram  at 
all,  say  over  L  adjacent  frequencies  (L=3  for  the  examples  in  this  paper),  the  distribution 
of  the  logarithm  of  the  smoothed  spectrum  will  be  close  to  a  normal  distribution.  This 
suggests  that,  as  an  approximation,  the  error  series  in  (21)  could  be  regarded  as  stationary 
normal  processes. 

BHT  (1962)  also  suggested  that,  in  the  case  of  a  single  series,  the  log  spectrum,  as  a 
function  of  frequency,  might  be  regarded  as  a  pseudo-time  series  and  that  the  cepstrum  be 
used  to  identify  the  periodicities  (time  delay  differences)  in  the  same  way  that  the  spectrum 
does  for  conventional  time  series.  In  accordance  with  previous  arguments,  it  is  clear  that 
|0(i^)|2  will  have  periods  equal  to  T/(tj  -  t*)  and  so  the  logarithm  of  the  periodogram  will 
tend  to  a  general  shape  having  peaks  at  (tj  -  Tk)/T,  j,  k  =  1, . . . ,  n.  Overall,  it  is  clear  that 
the  logarithm  of  the  periodogram  of  the  observed  series  in  (21)  will  still  be  nonstationary 
because  of  the  likelihood  that  the  component  log \Aj(u£)\2  will  vary  substantially  over 
frequency  and  in  a  manner  that  changes  for  different  recording  instruments. 
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Since  the  raw  series  is  also  expected  to  have  a  nonstationary  mean  component,  one 
could  detrend,  subtracting  a  reasonable  estimate  of  the  smooth  function  upon  which  the 
approximately  stationary  series  rides.  Any  relatively  smooth  detrending  will  work;  the 
method  used  here  was  to  fit  a  cubic  spline  with  a  single  knot  in  the  middle  of  the  record 
under  the  assumption  that  the  underlying  nonstationary  component  of  the  spectrum  is 
very  smooth  with  a  slight  inflection  point  about  half  way  through  the  record.  Moving 
average  smoothing,  called  liftering  by  BHT,  will  also  work  but  the  length  required  subjects 
the  series  to  serious  end  effects.  Once  the  detrending  has  been  done,  we  work  with  the 
resulting  series  and  do  not  make  further  use  of  the  functional  form  of  the  spline.  Here,  we 
envision  the  component  log  |A»|2  +  logP,(i/)  as  a  relatively  smooth  function  that  can 
be  approximated  by  a  cubic  spline  with  one  knot  as  in  (10).  We  estimate  the  parameters 
and  consider  the  detrended  log-periodogram  model 

\Y-(£)\2 

l0g  =  l0g  |e(l/<)|2  +  l0g  W Ml2-  (22) 

The  detrended  log-periodogram  is  expressed  in  terms  of  a  common  signal  log|0(i^)|2  and 
a  noise  term  log  |iVj(£)|2  which  differs  from  site  to  site  and  represents  local  site  effects. 

To  illustrate  some  of  the  above  ideas,  consider  again  the  simulated  ensemble  of  signals 
shown  in  the  top  panel  of  Figure  2.  The  five  resulting  series,  shown  in  Figure  2,  all 
contain  slightly  different  signals  at  the  same  set  of  delays  and  we  can  argue  that  the  model 
given  above  is  plausible  and  that  any  method  based  on  the  model  should  work  reasonably 
well.  A  typical  observed  spectrum  is  shown  in  the  bottom  panel  of  Figure  2  and  we  note 
that  it  appears  to  be  the  sum  of  a  smooth  trend  component,  log|AJ-(z')|2  +  log Ps(z/), 
shown  as  a  solid  line,  a  periodic  component,  corresponding  to  log|0(z/)|2  and  a  departure 
from  the  periodic  component  corresponding  to  log  |JVjz/)|2.  The  observed  log  spectrum 
is  nonstationary  because  of  the  smooth  component  and,  as  suggested  earlier,  it  would 
be  natural  to  detrend  the  log  spectrum  to  eliminate  the  effect  of  the  smooth  function 
log|Aj(i/)|2  +  log  Ps(u).  The  left  panel  of  Figure  5  shows  the  result  of  this  adjustment 
at  each  channel  as  dotted  lines  and  it  is  clear  that  one  might  regard  the  detrended  result 
as  the  sum  of  a  periodic  function  and  noise,  as  in  the  periodogram  model  given  above. 
Furthermore,  the  average  of  the  log  spectra  might  be  regarded  as  an  estimator  for  log  \9{v)  |2 
and  we  show  the  average  in  Figure  5  as  a  solid  line  superimposed  on  the  five  separate 
detrended  log  spectra.  The  average  seems  to  enhance  the  underlying  periodic  signal  and 
we  suggest  later  that  the  deviation  from  these  overall  means  are  basically  estimators  for 
the  log  spectra  of  the  noise  processes,  log  /”  (z/) .  The  result  of  detrending  with  a  cubic 
spline,  as  depicted  in  the  lower  panel  of  Figure  2,  has  the  appearance  of  a  stationary  signal 
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Figure  5:  Bank  of  Detrended  Log  Spectra  (Right)  and  the  Cepstral  F-Statistic  Show- 
ing  Peaks  at  8,  16,  23  and  30  points. 


with  several  periodic  components  corresponding  to  the  time  delay  differences  previously 
mentioned.  Here,  we  might  expect  periodicities  at  ri  -  rk  =  7,8,15,16,23,31,  with  po¬ 
tentially  the  strongest  components  corresponding  to  time  delays  of  7  or  8  points,  which 
appear  four  times,  time  delays  of  15  or  16  points,  which  appear  three  times,  a  time  delay 
of  23  points,  which  appears  twice  and  a  time  delay  of  31  points,  which  appears  once. 

This  example  and  model  (22)  above  motivates  an  approach  to  detecting  a  common 
signal  on  an  array  that  is  similar  to  that  given  in  Shumway  (1971),  where  the  discrete 
Fourier  transform,  say 


T/2-1 

Qj(d)  =  T-1/2  loS 
1=0 


WM 


exp{27r£&i/r} 


is  regarded  as  behaving  like  the  transform  of  a  stationary  process  in  the  delay  d  = 
0, 1,  •  •  • ,  T/4  at  each  channel  so  that  one  may  write  the  signal  plus  noise  model 


Qj(d)=S(d)  +  Vj(d),  (23) 

where  the  right  hand  side  contains  the  transforms  of  the  right  hand  side  of  (22).  Here, 
the  signal  transform  S  (d)  is  fixed  and  unknown  and  the  noise  Vj  (c/)  has  approximately  a 
complex  Gaussian  distribution  with  mean  0  and  variance  a2(d)  at  delay  d.  The  Fourier 
transform  of  such  a  process  will  be  in  the  form  of  a  sum  of  constants  multiplied  by  adjacent 
values  of  the  approximately  normally  distributed  log  periodogram  or  smoothed  spectra. 
The  linear  combination  will  be  approximately  normal,  since  the  components  are  approx¬ 
imately  normal,  and  will  tend  toward  normality  by  central  limit  arguments  even  when 
the  components  are  not  that  close  to  normality.  The  noises  are  assumed  uncorrelated  for 
different  j  and  have  identical  cepstra,  c2(d),  at  each  channel.  This  was  checked  for  several 
events  and  found  to  be  reasonable. 

Then,  motivated  by  the  classical  approach  to  detecting  a  signal  in  N  stationarily 
correlated  time  series,  we  note  that  testing  the  hypothesis  S(d)  =  0  leads  to  an  F-Statistic 
involving  the  total  stacked  cepstra  (SCT) 


SCTM  =  f;|C3,.(<i)|2  (24) 

3~  1 

and  the  spectrum  of  the  stacked  log  spectra  or  mean  stacked  cepstra  (SCM),  say 


SCM(d)  =  N\Q(d)\2, 


(25) 
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where 


(26) 


N 

Q(d)=N~1^QJ(d) 

j=l 

is  the  mean  Fourier  transform  of  the  array  log-spectra.  Figure  5  shows  the  original  de¬ 
trended  log-spectra  with  the  mean  of  the  array  log-spectra  superimposed  on  each  series  and 
we  note  that  the  common  periodicities  are  enhanced  in  the  stack.  An  important  quantity 
involved  in  the  optimal  detection  statistic  is  the  error  cepstrum  (SCE),  defined  as 

N 

SCE(d)  =  Y,\Qj(d)-Q(d)\2 

3- 1 

=  SCT{d)  -  SCM{d ),  (27) 


which  is  a  measure  of  the  extent  to  which  the  individual  channel  transforms  differ  from  the 
mean  transform.  It  can  be  interpreted  as  the  cepstral  noise  component.  The  F-Statistic 
resulting  from  the  signal  detection  hypothesis  is  given  by 


F2,2(N-i)(d)  =  (N-l) 


SCM  (d) 
SCE{d) 


(28) 


and  can  also  be  interpreted  as  a  cepstral  signal  to  noise  ratio.  The  subscripts  refer  to  an 
F  distribution  with  2  and  2  (AT  -  1)  degrees  of  freedom.  Note  that  the  total  cepstrum  (24) 
is  exactly  the  sum-stack  proposed  by  Alexander  et  al  (1995),  computed  by  adding  up  the 
separate  cepstra.  Alexander  et  al  have  also  considered  the  product-stack  which  does  not 
appear  to  have  any  identifiable  statistical  properties  and  we  do  not  analyze  it  here.  It  is 
clear  that  the  sum-stack  will  not  reflect  the  common  signal  components  as  well  as  either 
the  cepstral  component  due  to  the  signal  or  the  F-Statistic  (28). 

To  illustrate,  we  return  again  to  the  contrived  data  shown  in  Figure  2  which  contains 
signals  with  a  known  configuration  of  delays.  The  average  of  the  log  spectra,  shown  in 
the  left  panel  of  Figure  5,  is  an  approximately  unbiased  estimator  of  log  |#(i^)|2  in  (22); 
its  transform  estimates  S{d)  in  (23).  Figure  5  shows  the  F-Statistic  corresponding  to  the 
contrived  data  shown  in  Figure  2.  Note  the  strong  components  appearing  at  delays  of  8, 

15,  23,  30,  and  36  points  that  may  be  compared  with  the  known  true  delays  8,  15,  23,  and 
31  points.  Note  that  the  true  time  delays  would  imply  frequencies  of  the  form  7,  8,  15, 

16,  23,  and  31  points  respectively.  The  cepstral  F-Statistic,  shown  in  the  right  panel  of 
Figure  5,  provides  a  statistical  level  of  significance  for  the  various  peaks  and  we  note  that 
the  significant  peaks  are  8,  16,  23,  and  30  points  so  that  the  smallest  of  the  larger  peaks  at 
d=36  in  the  stacked  cepstrum  is  not  significant.  All  peaks  are  significant  at  a  false  alarm 
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rate  of  .001.  In  general,  since  there  are  often  a  large  number  of  delays  of  interest,  one 
should  insist  on  at  least  .01  as  a  level  of  significance. 

5.  Analysis  of  Kola  Mining  Explosions 

For  testing  on  actual  data,  we  focussed  on  the  deonvolution  and  cepstral  F  methodologies. 
The  parametric  seasonal  ARMA  models  of  Section  3  were  often  unstable  and  it  was  difficult 
to  obtain  reasonable  estimators  for  the  weights  9i,92,...,9n  for  the  simulated  data.  In 
addition,  the  estimation  procedure  seemed  to  be  unduly  sensitive  to  restricting  time  delays 
to  be  constant  multiples  of  some  underlying  delay  d.  In  contrast,  the  nonparametric 
deconvolution  and  cepstral  F  methods  concentrated  on  isolating  only  time  delay  differences 
on  the  theory  that  a  consistent  set  of  time  delay  differences  will  produce  a  strong  peak  in 
the  cepstrum  of  the  deconvolved  signal  or  in  the  cepstral  F-Statistic.  This  was  more  in 
line  with  the  idea  that  time  delays  induced  in  ripple-firing  may  be  somewhat  irregular.  For 
the  simulated  data,  both  frequency  domain  methods  isolated  the  predictable  time  delay 
differences. 

Data  used  for  testing  the  frequency  domain  methodology  were  from  mining  explosions 
in  the  Kola  Peninsula,  situated  in  the  Russian  Arctic.  Ten  mining  explosions  at  the  HD9 
quarry  were  observed  at  5  channels  of  the  ARCESS  Array  in  northern  Norway  (see  Der 
et  al,  1987).  Figure  1  shows  at  typical  mine  blast  (110),  sampled  at  .025  seconds,  and  we 
note  that  the  four  phases  Pn,Pg,Sn,  and  Lg  are  fairly  obvious.  This  is  typical  of  the  other 
events  (054,  066,  138,  147,  182,  219,  246,  282,  and  285)  and  it  was  fairly  easy  to  construct 
four  files  containing  the  separate  phases  for  each  of  the  10  events.  Figure  1  shows,  in  the 
bottom  panel,  the  five  series  used  for  the  Pn  phase  of  event  110. 

The  quantities  of  interest  will  be  the  cepstrum  of  the  deconvolved  signal  x(t),  given 
by  (11)  and  (12)  for  method  (i),  and  the  cepstral  F-Statistic  (28)  for  method  (iii).  Figures 
A1-A20  show  plots  of  the  two  statistics  for  each  of  the  10  events  (four  phases  per  plot). 
We  used  256  point  samples  for  all  phases  except  Lg,  which  had  512  points.  Ripple-firing 
may  be  indicated  best  by  a  consistent  pattern  of  strong  peaks,  with  delay  time  differences 
within  the  time  intervals  that  might  be  plausible  for  delays  or  their  aliases.  Sampling  is 
40  points  per  second  and  the  delays  are  given  in  points.  To  convert  to  milliseconds  (ms) 
multiply  by  25.  For  example,  Figure  Al  shows  peaks  at  4,  10-11,  16,  21,  and  26,  which 
might  imply  a  ripple  delay  of  5-6  points  or  of  125-150  ms.  Unfortunately,  Figure  A2, 
which  plots  the  F-Statistic  does  not  confirm  a  consistent  delay  in  that  range.  Note  that 
the  significance  values  for  the  F-Statistic  at  the  levels  .001,  .01,  and  .05  for  2  and  8  degrees 
of  freedom  are  18.5,  8.7,  and  4.5,  respectively. 


Table  1:  Delays  (in  pts.)  for  Deconvolved  Cepstra  and  Cepstral  F  Statistics 


Event 

Phase 

Pn 

P9 

Lg 

Deconvolutions 

054 

4,14,21,29 

5,12,21,26 

4,11,16,21,26 

5,15,28 

066 

2,10,23,29 

3,10,17,22 

8,14,20,28 

4,14 

110 

2,23 

7,17 

8,14,21,29 

3,9,18,25 

138 

3,9,14,23 

4,10,22,28 

4,11,31 

6,18,24,30 

147 

3,8,14,20,24,29 

6,17,21,26 

6,11,16,31 

3,7,11,15,20,29 

182 

4,14,21 

5,10 

7,19,26 

5,13,19,25 

219 

3,9,15,28 

2,8,11,15,22 

6,11,19 

3,7,11,15,23 

246 

3,8,13,18,27 

5,10,17,26 

3,7,11,20 

4,11,18,24 

282 

6,12,22 

2,15,20,25 

3,9,15,19,26 

2,11,16,26 

Cepstral  F 

.5in 

054 

4,21 

2,11,20 

2,12,21 

3 

066 

11,23 

11,15,22,30 

2,12 

110 

7,14,24 

,  8,25 

12 

2,9,21,25 

138 

2,8,18,31 

11,22,30 

2,12,19,26 

2,10,18,24 

147 

7,20 

7,31 

2,7,14,21,28 

2,7,15,24 

182 

4,21 

3,17,30 

8,13,20,28 

8,13,18,22,25 

219 

4,15,20,25,30 

6,15,26 

10,19,25 

8,11,16,23 

246 

12,23 

6,10,24,30 

7,13 

10,19,22 

282 

11,21,25,31 

14,22,29 

12,19,24,30 

3,11,16 

285 

3,11,17,22,26 

3,11,17,26 

5,13,18,27 

7,18 

We  went  through  all  the  plots  and  extracted  what  appeared  to  be  the  most  significant 
peaks  in  the  time  delays  for  both  the  cepstra  of  the  deconvolved  series  and  the  cepstral 
F-Statistic.  That  information  is  given  in  Table  1.  The  results  overall  are  somewhat  dis¬ 
appointing.  There  is  less  consistency  between  the  two  measures  than  one  would  have 
expected,  with  the  F-Statistic  generally  yielding  fewer  significant  peaks.  For  example,  Fig¬ 
ure  A8,  for  Event  138,  shows  evidence  of  weak  rippling  (.01  and  .05  levels)  in  the  Pn  and  Sn 
and  Lg  phases  but  not  in  Pg.  Figure  A7  from  the  same  event  shows  rippling  in  the  Pn  and 
Sn  phases  but  not  really  in  the  other  two  phases.  One  might  possibly  identify  consistent 
rippling  at  time  delays  from  6-10  points,  i.e.  150-250  ms  for  this  particular  event.  Other 
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events  are  roughly  similar,  with  the  evidence  of  delay  firing  often  indicated  weakly  in  the 
range  125-250  ms. 

It  is  clear  from  Table  1  that  all  phases  tend  to  show  evidence  of  delays  for  various 
events  but  that  there  are  no  clear  indications  that  the  same  patterns  for  all  phases  and  both 
statistics  hold  for  any  event  in  this  data  set.  Given  the  successes  enjoyed  by  both  methods 
on  simulated  data,  one  would  have  to  say  that  the  evidence  for  consistent  ripple  firing  over 
all  phases  in  the  Kola  mining  data  is  rather  weak.  In  addition,  there  is  not  a  single  phase 
that  seems  to  do  well  for  all  events.  Of  course,  a  larger  array,  say  of  30-40  elements,  could 
be  helpful  for  greater  enhancement  of  the  common  pattern  and  cancellation  of  patterns 
that  are  due  to  local  reverberations. 
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6.  Discussion 

In  this  project,  we  have  developed  and  tested  three  methods  for  detecting  ripple-firing  on 
arrays.  In  simulations,  it  is  clear  that  the  two  methods  based  on  the  cepstrum,  i.e.,  cepstral 
analysis  of  the  deconvolved  signal  and  the  cepstral  F-Statistic  offer  the  most  promise.  We 
prefer  the  cepstral  F-Statistic  for  the  following  reasons. 

1.  The  two  quantities  involved  in  the  cepstral  F  are  comparable  to  those  appearing  in 
the  cepstral  stack  proposed  by  Alexander  et  al  (1995)  for  the  purpose  of  estimating 
P-pP  reflection  delays.  Hence,  the  computed  statistic  is  a  function  of  quantities 
that  are  accepted  as  meaningful  in  the  seismic  literature. 

2.  The  statistic  has  a  convenient  interpretation  as  a  signal  to  noise  ratio  and  a 
known  statistical  distribution  under  the  hypothesis  of  no  delay  firing.  Hence,  we 
can  set  a  specific  level  for  false  alarms.  These  statistical  properties  distinguish  its 
performance  from  that  of  the  simple  stacked  spectra  of  Alexander  et  al  (1995). 

3.  Since  the  underlying  model  is  for  the  squared  periodogram,  it  is  not  critical  to 
fine  the  signals  up  on  the  array.  That  is,  the  form  is  related  to  incoherent  beam 
forming. 

It  would  be  remiss  not  to  mention  that  Alexander  et  al  (1995)  primarily  apply  their 
methods  to  the  problem  of  determining  source  depth.  This  is  functionally  related  to 
the  distance  between  the  primary  P  wave  and  the  single  reflection  pP.  This  problem  is 
somewhat  easier  since  there  is  conceptually  only  one  reflection  and  we  would  have  n  =  1 
in  (7),  i.e. 

\0(v)\2  =  1  +  $1  +  26i  cos(27ti >d) 

with  only  a  single  frequency  to  detect.  Hence,  it  would  be  interesting  to  apply  the  cepstral 
F-Statistic  to  the  depth  determination  problem  as  an  alternative  to  the  simple  sum  and 
product  stacks  employed  by  Alexander  et  al  (1995). 
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Cepstral  Analyses  of  Kola  Peninsula  Events 
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Figure  Al:  Cepstral  Analysis  of  Deconvolved  Array  Signal  for  Four  Phases  From 
Event  054  (Delay  in  points). 


Statistic  for  Four  Phases  From  Event  066  (Delay  in  points) 


PG110-F 


Figure  A7:  Cepstral  Analysis  of  Deconvolved  Array  Signal  for  Four  Phases  From 
Event  138  (Delay  in  points). 


PN138-F 


F-Statistic  for  Four  Phases  From  Event  138  (Delay  in  points). 
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stral  Analysis  of  Deconvolved  Array  Signal  for  Four  Phases  From 


Statistic  for  Four  Phases  From  Event  147  (Delay  in  points) 
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Statistic  for  Four  Phases  From  Event  246  (Delay  in  points). 


cepstrum  of  PN285  -  eepslrum  of  P6285 


Figure  A19:  Cepstral  Analysis  of  Deconvolved  Array  Signal  for  Four  Phases  From 
Event  285  (Delay  in  points). 


